Galaxy Liver Mets Kataoka et al_08052024 Clinical analysis

library(swimplot) library(coxphf) library(grid) library(gtable) library(readr) library(mosaic) library(dplyr) library(survival) library(survminer) library(gridtext) library(ggplot2) library(scales) library(officer) library(ggthemes) library(tidyverse) library(gtsummary) library(flextable) library(parameters) library(car) library(grid) library(ComplexHeatmap) library(readxl) library(janitor) library(rms) library(DT)

#Demographics Table

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data_subset <- circ_data %>%
  select(
    Age,
    Gender,
    ECOG,
    PrimSite,
    Mets.Type,
    Hepatectomy,
    NumLiverMetsGroup,
    SizeLiverMetsmmGroup,
    pT,
    pN,
    ACT,
    Oxaliplatin.History,
    Postop.Complication,
    BRAF.V600E,
    RAS,
    MSI,
    RFS.Event,
    OS.months) %>%
  mutate(
    Age = as.numeric(Age),
    Gender = factor(Gender, levels = c("Male", "Female")),
    ECOG = factor(ECOG, levels = c(0, 1)),
    PrimSite = factor(PrimSite, levels = c("Right-sided colon", "Left-sided colon")),
    Mets.Type = factor(Mets.Type, levels = c("Synchronous", "Metachronous")),
    Hepatectomy = factor(Hepatectomy, levels = c("Minor", "Major")),
    NumLiverMetsGroup = factor(NumLiverMetsGroup, levels = c("1", "≥2")),
    SizeLiverMetsmmGroup = factor(SizeLiverMetsmmGroup, levels = c("<50", "≥50")),
    pT = factor(pT, levels = c("T1-T2", "T3-T4")),
    pN = factor(pN, levels = c("N0", "N1-N2")),
    ACT = factor(ACT, levels = c("TRUE", "FALSE"), labels = c("Adjuvant Chemotherapy", "Observation")),
    Oxaliplatin.History = factor(Oxaliplatin.History, levels = c("No", "Yes")),
    Postop.Complication = factor(Postop.Complication, levels = c("No", "Yes")),
    BRAF.V600E = factor(BRAF.V600E, levels = c("WT", "MUT"), labels = c("BRAF wt", "BRAF V600E")),
    RAS = factor(RAS, levels = c("WT", "MUT"), labels = c("RAS wt", "RAS mut")),
    MSI = factor(MSI, levels = c("MSS", "MSI-High")),
    RFS.Event = factor(RFS.Event, levels = c("TRUE", "FALSE"), labels = c("Recurrence", "No Recurrence")),
    OS.months = as.numeric(OS.months))
table1 <- circ_data_subset %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  bold_labels()
table1

Characteristic	N = 1901
Age	68 (34 - 85)
Gender
Male	118 (62%)
Female	72 (38%)
ECOG
0	186 (98%)
1	4 (2.1%)
PrimSite
Right-sided colon	50 (26%)
Left-sided colon	140 (74%)
Mets.Type
Synchronous	64 (34%)
Metachronous	126 (66%)
Hepatectomy
Minor	159 (84%)
Major	31 (16%)
NumLiverMetsGroup
1	121 (64%)
≥2	69 (36%)
SizeLiverMetsmmGroup
<50	180 (95%)
≥50	10 (5.3%)
pT
T1-T2	23 (12%)
T3-T4	167 (88%)
pN
N0	77 (41%)
N1-N2	113 (59%)
ACT
Adjuvant Chemotherapy	48 (25%)
Observation	142 (75%)
Oxaliplatin.History	50 (26%)
Postop.Complication	23 (12%)
BRAF.V600E
BRAF wt	188 (99%)
BRAF V600E	2 (1.1%)
RAS
RAS wt	106 (56%)
RAS mut	84 (44%)
MSI
MSS	188 (99%)
MSI-High	2 (1.1%)
RFS.Event
Recurrence	95 (50%)
No Recurrence	95 (50%)
OS.months	24 (1 - 48)
1 Median (Range); n (%)

fit1 <- as_flex_table(
  table1,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)

Warning: The `strip_md_bold` argument of `as_flex_table()` is deprecated as of gtsummary 1.6.0.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

fit1

Characteristic	N = 1901
Age	68 (34 - 85)
Gender
Male	118 (62%)
Female	72 (38%)
ECOG
0	186 (98%)
1	4 (2.1%)
PrimSite
Right-sided colon	50 (26%)
Left-sided colon	140 (74%)
Mets.Type
Synchronous	64 (34%)
Metachronous	126 (66%)
Hepatectomy
Minor	159 (84%)
Major	31 (16%)
NumLiverMetsGroup
1	121 (64%)
≥2	69 (36%)
SizeLiverMetsmmGroup
<50	180 (95%)
≥50	10 (5.3%)
pT
T1-T2	23 (12%)
T3-T4	167 (88%)
pN
N0	77 (41%)
N1-N2	113 (59%)
ACT
Adjuvant Chemotherapy	48 (25%)
Observation	142 (75%)
Oxaliplatin.History	50 (26%)
Postop.Complication	23 (12%)
BRAF.V600E
BRAF wt	188 (99%)
BRAF V600E	2 (1.1%)
RAS
RAS wt	106 (56%)
RAS mut	84 (44%)
MSI
MSS	188 (99%)
MSI-High	2 (1.1%)
RFS.Event
Recurrence	95 (50%)
No Recurrence	95 (50%)
OS.months	24 (1 - 48)
1Median (Range); n (%)

save_as_docx(fit1, path= "~/Downloads/table1.docx")

#Demographics Table by MRD ctDNA Status

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data_subset1 <- circ_data %>%
  select(
    Age,
    Gender,
    ECOG,
    PrimSite,
    Mets.Type,
    Hepatectomy,
    NumLiverMetsGroup,
    SizeLiverMetsmmGroup,
    pT,
    pN,
    ACT,
    Oxaliplatin.History,
    Postop.Complication,
    BRAF.V600E,
    RAS,
    MSI) %>%
  mutate(
    Age = as.numeric(Age),
    Gender = factor(Gender, levels = c("Male", "Female")),
    ECOG = factor(ECOG, levels = c(0, 1)),
    PrimSite = factor(PrimSite, levels = c("Right-sided colon", "Left-sided colon")),
    Mets.Type = factor(Mets.Type, levels = c("Synchronous", "Metachronous")),
    Hepatectomy = factor(Hepatectomy, levels = c("Minor", "Major")),
    NumLiverMetsGroup = factor(NumLiverMetsGroup, levels = c("1", "≥2")),
    SizeLiverMetsmmGroup = factor(SizeLiverMetsmmGroup, levels = c("<50", "≥50")),
    pT = factor(pT, levels = c("T1-T2", "T3-T4")),
    pN = factor(pN, levels = c("N0", "N1-N2")),
    ACT = factor(ACT, levels = c("TRUE", "FALSE"), labels = c("Adjuvant Chemotherapy", "Observation")),
    Oxaliplatin.History = factor(Oxaliplatin.History, levels = c("No", "Yes")),
    Postop.Complication = factor(Postop.Complication, levels = c("No", "Yes")),
    BRAF.V600E = factor(BRAF.V600E, levels = c("WT", "MUT"), labels = c("BRAF wt", "BRAF V600E")),
    RAS = factor(RAS, levels = c("WT", "MUT"), labels = c("RAS wt", "RAS mut")),
    MSI = factor(MSI, levels = c("MSS", "MSI-High")))

circ_data1 <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data1 <- circ_data1[circ_data1$LiverMets=="TRUE",]

circ_data_subset2 <- circ_data %>%
  select(
    Age,
    Gender,
    ECOG,
    PrimSite,
    Mets.Type,
    Hepatectomy,
    NumLiverMetsGroup,
    SizeLiverMetsmmGroup,
    pT,
    pN,
    ACT,
    Oxaliplatin.History,
    Postop.Complication,
    BRAF.V600E,
    RAS,
    MSI,
    ctDNA.MRD) %>%
  mutate(
    Age = as.numeric(Age),
    Gender = factor(Gender, levels = c("Male", "Female")),
    ECOG = factor(ECOG, levels = c(0, 1)),
    PrimSite = factor(PrimSite, levels = c("Right-sided colon", "Left-sided colon")),
    Mets.Type = factor(Mets.Type, levels = c("Synchronous", "Metachronous")),
    Hepatectomy = factor(Hepatectomy, levels = c("Minor", "Major")),
    NumLiverMetsGroup = factor(NumLiverMetsGroup, levels = c("1", "≥2")),
    SizeLiverMetsmmGroup = factor(SizeLiverMetsmmGroup, levels = c("<50", "≥50")),
    pT = factor(pT, levels = c("T1-T2", "T3-T4")),
    pN = factor(pN, levels = c("N0", "N1-N2")),
    ACT = factor(ACT, levels = c("TRUE", "FALSE"), labels = c("Adjuvant Chemotherapy", "Observation")),
    Oxaliplatin.History = factor(Oxaliplatin.History, levels = c("No", "Yes")),
    Postop.Complication = factor(Postop.Complication, levels = c("No", "Yes")),
    BRAF.V600E = factor(BRAF.V600E, levels = c("WT", "MUT"), labels = c("BRAF wt", "BRAF V600E")),
    RAS = factor(RAS, levels = c("WT", "MUT"), labels = c("RAS wt", "RAS mut")),
    MSI = factor(MSI, levels = c("MSS", "MSI-High")),
    ctDNA.MRD = factor(ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE")))
Overall <- circ_data_subset1 %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  bold_labels()
Overall

Characteristic	N = 1901
Age	68 (34 - 85)
Gender
Male	118 (62%)
Female	72 (38%)
ECOG
0	186 (98%)
1	4 (2.1%)
PrimSite
Right-sided colon	50 (26%)
Left-sided colon	140 (74%)
Mets.Type
Synchronous	64 (34%)
Metachronous	126 (66%)
Hepatectomy
Minor	159 (84%)
Major	31 (16%)
NumLiverMetsGroup
1	121 (64%)
≥2	69 (36%)
SizeLiverMetsmmGroup
<50	180 (95%)
≥50	10 (5.3%)
pT
T1-T2	23 (12%)
T3-T4	167 (88%)
pN
N0	77 (41%)
N1-N2	113 (59%)
ACT
Adjuvant Chemotherapy	48 (25%)
Observation	142 (75%)
Oxaliplatin.History	50 (26%)
Postop.Complication	23 (12%)
BRAF.V600E
BRAF wt	188 (99%)
BRAF V600E	2 (1.1%)
RAS
RAS wt	106 (56%)
RAS mut	84 (44%)
MSI
MSS	188 (99%)
MSI-High	2 (1.1%)
1 Median (Range); n (%)

ByctDNA_MRD <- circ_data_subset2 %>%
  tbl_summary(
    by = ctDNA.MRD, # Subgroup by ctDNA.MRD
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  add_p(test = all_categorical() ~ "fisher.test", pvalue_fun = ~style_pvalue(p.adjust(.x, method = "bonferroni"))) %>%
  bold_labels()
ByctDNA_MRD

Characteristic	NEGATIVE, N = 1291	POSITIVE, N = 611	p-value2
Age	69 (37 - 85)	68 (34 - 84)	>0.9
Gender			>0.9
Male	83 (64%)	35 (57%)
Female	46 (36%)	26 (43%)
ECOG			>0.9
0	127 (98%)	59 (97%)
1	2 (1.6%)	2 (3.3%)
PrimSite			0.6
Right-sided colon	40 (31%)	10 (16%)
Left-sided colon	89 (69%)	51 (84%)
Mets.Type			0.015
Synchronous	33 (26%)	31 (51%)
Metachronous	96 (74%)	30 (49%)
Hepatectomy			>0.9
Minor	112 (87%)	47 (77%)
Major	17 (13%)	14 (23%)
NumLiverMetsGroup			0.010
1	93 (72%)	28 (46%)
≥2	36 (28%)	33 (54%)
SizeLiverMetsmmGroup			>0.9
<50	123 (95%)	57 (93%)
≥50	6 (4.7%)	4 (6.6%)
pT			>0.9
T1-T2	16 (12%)	7 (11%)
T3-T4	113 (88%)	54 (89%)
pN			0.3
N0	60 (47%)	17 (28%)
N1-N2	69 (53%)	44 (72%)
ACT			>0.9
Adjuvant Chemotherapy	33 (26%)	15 (25%)
Observation	96 (74%)	46 (75%)
Oxaliplatin.History	39 (30%)	11 (18%)	>0.9
Postop.Complication	10 (7.8%)	13 (21%)	0.2
BRAF.V600E			>0.9
BRAF wt	127 (98%)	61 (100%)
BRAF V600E	2 (1.6%)	0 (0%)
RAS			>0.9
RAS wt	74 (57%)	32 (52%)
RAS mut	55 (43%)	29 (48%)
MSI			>0.9
MSS	128 (99%)	60 (98%)
MSI-High	1 (0.8%)	1 (1.6%)
1 Median (Range); n (%)
2 Wilcoxon rank sum test; Fisher’s exact test

merged_table <- tbl_merge(tbls=list(Overall, ByctDNA_MRD))
merged_table

Characteristic	Table 1	Table 2
	N = 1901	NEGATIVE, N = 1291	POSITIVE, N = 611	p-value2
Age	68 (34 - 85)	69 (37 - 85)	68 (34 - 84)	>0.9
Gender				>0.9
Male	118 (62%)	83 (64%)	35 (57%)
Female	72 (38%)	46 (36%)	26 (43%)
ECOG				>0.9
0	186 (98%)	127 (98%)	59 (97%)
1	4 (2.1%)	2 (1.6%)	2 (3.3%)
PrimSite				0.6
Right-sided colon	50 (26%)	40 (31%)	10 (16%)
Left-sided colon	140 (74%)	89 (69%)	51 (84%)
Mets.Type				0.015
Synchronous	64 (34%)	33 (26%)	31 (51%)
Metachronous	126 (66%)	96 (74%)	30 (49%)
Hepatectomy				>0.9
Minor	159 (84%)	112 (87%)	47 (77%)
Major	31 (16%)	17 (13%)	14 (23%)
NumLiverMetsGroup				0.010
1	121 (64%)	93 (72%)	28 (46%)
≥2	69 (36%)	36 (28%)	33 (54%)
SizeLiverMetsmmGroup				>0.9
<50	180 (95%)	123 (95%)	57 (93%)
≥50	10 (5.3%)	6 (4.7%)	4 (6.6%)
pT				>0.9
T1-T2	23 (12%)	16 (12%)	7 (11%)
T3-T4	167 (88%)	113 (88%)	54 (89%)
pN				0.3
N0	77 (41%)	60 (47%)	17 (28%)
N1-N2	113 (59%)	69 (53%)	44 (72%)
ACT				>0.9
Adjuvant Chemotherapy	48 (25%)	33 (26%)	15 (25%)
Observation	142 (75%)	96 (74%)	46 (75%)
Oxaliplatin.History	50 (26%)	39 (30%)	11 (18%)	>0.9
Postop.Complication	23 (12%)	10 (7.8%)	13 (21%)	0.2
BRAF.V600E				>0.9
BRAF wt	188 (99%)	127 (98%)	61 (100%)
BRAF V600E	2 (1.1%)	2 (1.6%)	0 (0%)
RAS				>0.9
RAS wt	106 (56%)	74 (57%)	32 (52%)
RAS mut	84 (44%)	55 (43%)	29 (48%)
MSI				>0.9
MSS	188 (99%)	128 (99%)	60 (98%)
MSI-High	2 (1.1%)	1 (0.8%)	1 (1.6%)
1 Median (Range); n (%)
2 Wilcoxon rank sum test; Fisher’s exact test

fit1 <- as_flex_table(
  merged_table,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)
fit1

	Table 1	Table 2
Characteristic	N = 1901	NEGATIVE, N = 1291	POSITIVE, N = 611	p-value2
Age	68 (34 - 85)	69 (37 - 85)	68 (34 - 84)	>0.9
Gender				>0.9
Male	118 (62%)	83 (64%)	35 (57%)
Female	72 (38%)	46 (36%)	26 (43%)
ECOG				>0.9
0	186 (98%)	127 (98%)	59 (97%)
1	4 (2.1%)	2 (1.6%)	2 (3.3%)
PrimSite				0.6
Right-sided colon	50 (26%)	40 (31%)	10 (16%)
Left-sided colon	140 (74%)	89 (69%)	51 (84%)
Mets.Type				0.015
Synchronous	64 (34%)	33 (26%)	31 (51%)
Metachronous	126 (66%)	96 (74%)	30 (49%)
Hepatectomy				>0.9
Minor	159 (84%)	112 (87%)	47 (77%)
Major	31 (16%)	17 (13%)	14 (23%)
NumLiverMetsGroup				0.010
1	121 (64%)	93 (72%)	28 (46%)
≥2	69 (36%)	36 (28%)	33 (54%)
SizeLiverMetsmmGroup				>0.9
<50	180 (95%)	123 (95%)	57 (93%)
≥50	10 (5.3%)	6 (4.7%)	4 (6.6%)
pT				>0.9
T1-T2	23 (12%)	16 (12%)	7 (11%)
T3-T4	167 (88%)	113 (88%)	54 (89%)
pN				0.3
N0	77 (41%)	60 (47%)	17 (28%)
N1-N2	113 (59%)	69 (53%)	44 (72%)
ACT				>0.9
Adjuvant Chemotherapy	48 (25%)	33 (26%)	15 (25%)
Observation	142 (75%)	96 (74%)	46 (75%)
Oxaliplatin.History	50 (26%)	39 (30%)	11 (18%)	>0.9
Postop.Complication	23 (12%)	10 (7.8%)	13 (21%)	0.2
BRAF.V600E				>0.9
BRAF wt	188 (99%)	127 (98%)	61 (100%)
BRAF V600E	2 (1.1%)	2 (1.6%)	0 (0%)
RAS				>0.9
RAS wt	106 (56%)	74 (57%)	32 (52%)
RAS mut	84 (44%)	55 (43%)	29 (48%)
MSI				>0.9
MSS	188 (99%)	128 (99%)	60 (98%)
MSI-High	2 (1.1%)	1 (0.8%)	1 (1.6%)
1Median (Range); n (%)
2Wilcoxon rank sum test; Fisher's exact test

save_as_docx(fit1, path= "~/Downloads/merged_table.docx")

#ctDNA Detection Rates by Window and Stages

#ctDNA at Baseline
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Baseline!="",]
circ_data$ctDNA.Baseline <- factor(circ_data$ctDNA.Baseline, levels=c("NEGATIVE","POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Baseline == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Baseline, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Baseline == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA at MRD Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.MRD == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.MRD, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.MRD == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA at Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
circ_data <- subset(circ_data, ctDNA.Surveillance %in% c("NEGATIVE", "POSITIVE"))
circ_data$Stage <- factor(circ_data$Stage, levels=c("IV"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance == "POSITIVE", by=list(circ_data$Stage), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance, by=list(circ_data$Stage), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Surveillance == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA Detection Rates by Window and Synchronicity

#ctDNA at Baseline
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Baseline!="",]
circ_data$ctDNA.Baseline <- factor(circ_data$ctDNA.Baseline, levels=c("NEGATIVE","POSITIVE"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Baseline == "POSITIVE", by=list(circ_data$Mets.Type), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Baseline, by=list(circ_data$Mets.Type), FUN=length)
combined_data <- data.frame(
  Stage = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Baseline == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Stage = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA at MRD Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.MRD == "POSITIVE", by=list(circ_data$Mets.Type), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.MRD, by=list(circ_data$Mets.Type), FUN=length)
combined_data <- data.frame(
  Mets.Type = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.MRD == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Mets.Type = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#ctDNA at Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
circ_data <- subset(circ_data, ctDNA.Surveillance %in% c("NEGATIVE", "POSITIVE"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
positive_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance == "POSITIVE", by=list(circ_data$Mets.Type), FUN=sum)
total_counts_by_stage <- aggregate(circ_data$ctDNA.Surveillance, by=list(circ_data$Mets.Type), FUN=length)
combined_data <- data.frame(
  Mets.Type = total_counts_by_stage$Group.1,
  Total_Count = total_counts_by_stage$x,
  Positive_Count = positive_counts_by_stage$x,
  Rate = (positive_counts_by_stage$x / total_counts_by_stage$x) * 100  # Convert to percentage
)
combined_data$Rate <- sprintf("%.2f%%", combined_data$Rate)
overall_total_count <- nrow(circ_data)
overall_positive_count <- nrow(circ_data[circ_data$ctDNA.Surveillance == "POSITIVE",])
overall_positivity_rate <- (overall_positive_count / overall_total_count) * 100  # Convert to percentage
overall_row <- data.frame(
  Mets.Type = "Overall",
  Total_Count = overall_total_count,
  Positive_Count = overall_positive_count,
  Rate = sprintf("%.2f%%", overall_positivity_rate)
)
combined_data <- rbind(combined_data, overall_row)
print(combined_data)

#Barplot with Detection Rates at the MRD and Surveillance Windows by Synchronicity

#Detection rate at the Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
contingency_table <- table(circ_data$Mets.Type, circ_data$ctDNA.MRD)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 10.707, df = 1, p-value = 0.001067

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.0009273
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.507768 5.977774
sample estimates:
odds ratio 
  2.987004 

print(contingency_table)

              
               Negative Positive
  Metachronous       96       30
  Synchronous        33       31

p_values <- c(chi_square_test$p.value, fisher_exact_test$p.value)
p_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_adjusted) <- c("Chi-Square Test", "Fisher's Exact Test")
print(p_adjusted)

    Chi-Square Test Fisher's Exact Test 
        0.002134653         0.001854561 

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Patients with CLM - ctDNA MRD Window", 
       x = "Synchronicity", 
       y = "Patients (%)", 
       fill = "ctDNA MRD Window",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value), "; Bonferroni corrected p-value: ", format.pval(p_adjusted["Chi-Square Test"]))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#Detection rate at the Surveillance Window
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
contingency_table <- table(circ_data$Mets.Type, circ_data$ctDNA.Surveillance)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 4.4449, df = 1, p-value = 0.03501

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.03232
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.048856 4.688132
sample estimates:
odds ratio 
  2.210114 

print(contingency_table)

              
               Negative Positive
  Metachronous       72       30
  Synchronous        27       25

p_values <- c(chi_square_test$p.value, fisher_exact_test$p.value)
p_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_adjusted) <- c("Chi-Square Test", "Fisher's Exact Test")
print(p_adjusted)

    Chi-Square Test Fisher's Exact Test 
         0.07001175          0.06464745 

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Patients with CLM - ctDNA Surveillance Window", 
       x = "Synchronicity", 
       y = "Patients (%)", 
       fill = "ctDNA Surveillance Window",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value), "; Bonferroni corrected p-value: ", format.pval(p_adjusted["Chi-Square Test"]))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#Barplot with Detection Rates at pre-op and ctDNA MRD by pre-OP CEA

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data$ctDNA.Baseline <- factor(circ_data$ctDNA.Baseline, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$CEA.Baseline <- factor(circ_data$CEA.Baseline, levels=c("NEGATIVE","POSITIVE"), labels = c("Normal", "Elevated"))
contingency_table <- table(circ_data$CEA.Baseline, circ_data$ctDNA.Baseline)
chi_square_test <- chisq.test(contingency_table)

Warning in stats::chisq.test(x, y, ...) :
  Chi-squared approximation may be incorrect

print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 0.038105, df = 1, p-value = 0.8452

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 0.5858
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
   0.1293152 150.9263898
sample estimates:
odds ratio 
  2.524321 

print(contingency_table)

          
           Negative Positive
  Normal          2       82
  Elevated        1      104

p_values <- c(chi_square_test$p.value, fisher_exact_test$p.value)
p_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_adjusted) <- c("Chi-Square Test", "Fisher's Exact Test")
print(p_adjusted)

    Chi-Square Test Fisher's Exact Test 
                  1                   1 

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Pre-op CEA & Pre-op ctDNA", 
       x = "Pre-op CEA", 
       y = "Patients (%)", 
       fill = "ctDNA pre-op",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value), "; Bonferroni corrected p-value: ", format.pval(p_adjusted["Chi-Square Test"]))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$CEA.Baseline <- factor(circ_data$CEA.Baseline, levels=c("NEGATIVE","POSITIVE"), labels = c("Normal", "Elevated"))
contingency_table <- table(circ_data$CEA.Baseline, circ_data$ctDNA.MRD)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 2.8368e-30, df = 1, p-value = 1

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 1
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.5168345 1.9328504
sample estimates:
odds ratio 
 0.9969296 

print(contingency_table)

          
           Negative Positive
  Normal         57       27
  Elevated       72       34

p_values <- c(chi_square_test$p.value, fisher_exact_test$p.value)
p_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_adjusted) <- c("Chi-Square Test", "Fisher's Exact Test")
print(p_adjusted)

    Chi-Square Test Fisher's Exact Test 
                  1                   1 

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Pre-op CEA & ctDNA MRD", 
       x = "Pre-op CEA", 
       y = "Patients (%)", 
       fill = "ctDNA at MRD",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value), "; Bonferroni corrected p-value: ", format.pval(p_adjusted["Chi-Square Test"]))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#DFS by ctDNA at the MRD Window - Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                     n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 129     45     NA      32      NA
ctDNA.MRD=POSITIVE  60     53   2.56       2     4.6

event_summary <- circ_data %>%
  group_by(ctDNA.MRD) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA MRD window", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     41      42    0.645  0.0456        0.548        0.726
   30     29       2    0.613  0.0486        0.511        0.700
   36     15       1    0.592  0.0513        0.484        0.684

                ctDNA.MRD=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24      4      53    0.108  0.0409       0.0445        0.203
   30      3       0    0.108  0.0409       0.0445        0.203
   36      2       0    0.108  0.0409       0.0445        0.203

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 189, number of events= 98 

                    coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRDPOSITIVE 1.7470    5.7371   0.2092 8.352   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     5.737     0.1743     3.808     8.645

Concordance= 0.701  (se = 0.021 )
Likelihood ratio test= 65.87  on 1 df,   p=5e-16
Wald test            = 69.75  on 1 df,   p=<2e-16
Score (logrank) test = 86.49  on 1 df,   p=<2e-16

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 5.74 (3.81-8.64); p = 0"

#Fisher test for DFS percentages at 24, 30, and 36 months
dfs_times <- c(24, 30, 36)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$DFS.MRD.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$DFS.MRD.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.MRD.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.MRD.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
        2.248794e-13         8.823505e-13         1.412511e-12 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", dfs_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months Corrected p-value at 30 months Corrected p-value at 36 months 
                  6.746383e-13                   2.647051e-12                   4.237533e-12 

#OS by ctDNA at the MRD Window - Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$OS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                     n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 129      5     NA      NA      NA
ctDNA.MRD=POSITIVE  61     11   43.4      NA      NA

event_summary <- circ_data %>%
  group_by(ctDNA.MRD) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA MRD window | All stages", ylab= "Overall Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     71       2    0.984  0.0111        0.938        0.996
   30     48       0    0.984  0.0111        0.938        0.996
   36     26       2    0.943  0.0306        0.841        0.980

                ctDNA.MRD=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     23       7    0.836  0.0601        0.674        0.921
   30     12       3    0.717  0.0819        0.521        0.844
   36      6       0    0.717  0.0819        0.521        0.844

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 190, number of events= 16 

                    coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.MRDPOSITIVE 1.8630    6.4429   0.5441 3.424 0.000618 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     6.443     0.1552     2.218     18.72

Concordance= 0.726  (se = 0.064 )
Likelihood ratio test= 12.93  on 1 df,   p=3e-04
Wald test            = 11.72  on 1 df,   p=6e-04
Score (logrank) test = 15.32  on 1 df,   p=9e-05

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 4))
print(label_text)

[1] "HR = 6.44 (2.22-18.72); p = 6e-04"

#Fisher test for OS percentages at 24, 30, and 36 months
os_times <- c(24, 30, 36)
p_values <- sapply(os_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", os_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
        0.0054510839         0.0002468767         0.0020570349 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", os_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months Corrected p-value at 30 months Corrected p-value at 36 months 
                  0.0163532517                   0.0007406301                   0.0061711046 

#Multivariate cox regression at MRD Window for DFS - Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$CEA.Baseline <- factor(circ_data$CEA.Baseline, levels=c("NEGATIVE","POSITIVE"), labels = c("Normal", "Elevated"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.MRD + CEA.Baseline + Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression at MRD Window for DFS - ctDNA MRD Positive Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$CEA.Baseline <- factor(circ_data$CEA.Baseline, levels=c("NEGATIVE","POSITIVE"), labels = c("Normal", "Elevated"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ CEA.Baseline + Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS - ctDNA MRD Positive", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#Multivariate cox regression at MRD Window for DFS - ctDNA MRD Negative Landmark MRD timepoint

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$DFS.MRD.months>=0,]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_datadf <- as.data.frame(circ_data)

circ_data$CEA.Baseline <- factor(circ_data$CEA.Baseline, levels=c("NEGATIVE","POSITIVE"), labels = c("Normal", "Elevated"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event)
cox_fit <- coxph(surv_object ~ CEA.Baseline + Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS - ctDNA MRD Negative", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#DFS by ACT treatment in MRD negative

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 92     33     NA    25.6      NA
ACT=TRUE  33      8     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      25.0000      31.0000       0.6222       0.0566       0.5011       0.7218 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      15.0000       8.0000       0.7231       0.0859       0.5143       0.8538 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 125, number of events= 41 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.5575    1.7463   0.3948 1.412    0.158

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     1.746     0.5727    0.8054     3.786

Concordance= 0.553  (se = 0.032 )
Likelihood ratio test= 2.22  on 1 df,   p=0.1
Wald test            = 1.99  on 1 df,   p=0.2
Score (logrank) test = 2.05  on 1 df,   p=0.2

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.75 (0.81-3.79); p = 0.158"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
           0.3844277 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", dfs_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months 
                     0.3844277 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, 
    data = circ_data)

  n= 125, number of events= 41 

                           coef exp(coef) se(coef)      z Pr(>|z|)
ACTTRUE                 -0.3844    0.6809   0.4300 -0.894    0.371
GenderMale               0.3311    1.3925   0.3938  0.841    0.400
Age.Group>70             0.2674    1.3066   0.3279  0.815    0.415
NumLiverMetsGroup≥2      0.2951    1.3432   0.3658  0.807    0.420
SizeLiverMetsmmGroup≥50  0.4474    1.5642   0.6244  0.716    0.474
RASMutant                0.3107    1.3644   0.3582  0.867    0.386
Mets.TypeSynchronous    -0.1946    0.8232   0.4471 -0.435    0.663
Oxaliplatin.HistoryYes   0.2572    1.2934   0.3778  0.681    0.496

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.6809     1.4687    0.2931     1.581
GenderMale                 1.3925     0.7181    0.6436     3.013
Age.Group>70               1.3066     0.7654    0.6871     2.485
NumLiverMetsGroup≥2        1.3432     0.7445    0.6559     2.751
SizeLiverMetsmmGroup≥50    1.5642     0.6393    0.4600     5.319
RASMutant                  1.3644     0.7329    0.6761     2.753
Mets.TypeSynchronous       0.8232     1.2148    0.3427     1.977
Oxaliplatin.HistoryYes     1.2934     0.7732    0.6168     2.712

Concordance= 0.638  (se = 0.042 )
Likelihood ratio test= 6.52  on 8 df,   p=0.6
Wald test            = 6.23  on 8 df,   p=0.6
Score (logrank) test = 6.37  on 8 df,   p=0.6

#DFS by ACT treatment in MRD positive

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 31     29   1.45    0.76    3.32
ACT=TRUE  15     10  12.92    4.11      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Positive ACT vs Observation", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(12, 24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     12.0000       3.0000      27.0000       0.1290       0.0602       0.0407       0.2698 

                ACT=TRUE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   12      9       6    0.600   0.126        0.318        0.797
   24      4       4    0.333   0.122        0.122        0.564

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 46, number of events= 39 

           coef exp(coef) se(coef)     z Pr(>|z|)    
ACTFALSE 1.3838    3.9899   0.3842 3.602 0.000316 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE      3.99     0.2506     1.879     8.472

Concordance= 0.675  (se = 0.03 )
Likelihood ratio test= 15  on 1 df,   p=1e-04
Wald test            = 12.97  on 1 df,   p=3e-04
Score (logrank) test = 14.61  on 1 df,   p=1e-04

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 3.99 (1.88-8.47); p = 0"

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Mets.Type + Oxaliplatin.History, 
    data = circ_data)

  n= 46, number of events= 39 

                            coef exp(coef) se(coef)      z Pr(>|z|)    
ACTTRUE                 -2.68676   0.06810  0.67360 -3.989 6.65e-05 ***
GenderMale              -1.04107   0.35308  0.51890 -2.006  0.04483 *  
Age.Group>70            -0.01813   0.98203  0.39282 -0.046  0.96318    
NumLiverMetsGroup≥2     -0.94250   0.38965  0.43212 -2.181  0.02918 *  
SizeLiverMetsmmGroup≥50  1.96295   7.12027  0.91333  2.149  0.03162 *  
RASMutant                1.01046   2.74687  0.36487  2.769  0.00562 ** 
Mets.TypeSynchronous     0.62120   1.86116  0.44165  1.407  0.15956    
Oxaliplatin.HistoryYes   0.81261   2.25379  0.57133  1.422  0.15493    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.0681    14.6840   0.01819    0.2550
GenderMale                 0.3531     2.8322   0.12770    0.9762
Age.Group>70               0.9820     1.0183   0.45473    2.1208
NumLiverMetsGroup≥2        0.3897     2.5664   0.16705    0.9089
SizeLiverMetsmmGroup≥50    7.1203     0.1404   1.18869   42.6504
RASMutant                  2.7469     0.3641   1.34357    5.6159
Mets.TypeSynchronous       1.8612     0.5373   0.78316    4.4230
Oxaliplatin.HistoryYes     2.2538     0.4437   0.73552    6.9061

Concordance= 0.785  (se = 0.033 )
Likelihood ratio test= 35.06  on 8 df,   p=3e-05
Wald test            = 28.56  on 8 df,   p=4e-04
Score (logrank) test = 35.96  on 8 df,   p=2e-05

#DFS by ACT treatment in MRD negative - Synchronous metastatic site

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Synchronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 16      5     NA    9.01      NA
ACT=TRUE  17      3     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | Synchronous Met Site", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        5.000        5.000        0.661        0.125        0.364        0.844 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        9.000        3.000        0.776        0.116        0.447        0.924 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 33, number of events= 8 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.8259    2.2838   0.7334 1.126     0.26

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     2.284     0.4379    0.5425     9.615

Concordance= 0.634  (se = 0.082 )
Likelihood ratio test= 1.32  on 1 df,   p=0.3
Wald test            = 1.27  on 1 df,   p=0.3
Score (logrank) test = 1.34  on 1 df,   p=0.2

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 2.28 (0.54-9.62); p = 0.26"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
           0.4384333 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", dfs_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months 
                     0.4384333 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Synchronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 33, number of events= 8 

                           coef exp(coef) se(coef)      z Pr(>|z|)  
ACTTRUE                 -1.5633    0.2094   0.8370 -1.868   0.0618 .
GenderMale              -1.1347    0.3215   0.9708 -1.169   0.2425  
Age.Group>70            -0.6898    0.5017   0.8681 -0.795   0.4268  
NumLiverMetsGroup≥2      1.0846    2.9583   0.9082  1.194   0.2324  
SizeLiverMetsmmGroup≥50  1.7630    5.8302   1.3825  1.275   0.2022  
RASMutant                1.5951    4.9288   1.0234  1.559   0.1191  
Oxaliplatin.HistoryYes       NA        NA   0.0000     NA       NA  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.2094     4.7744   0.04061     1.080
GenderMale                 0.3215     3.1102   0.04796     2.155
Age.Group>70               0.5017     1.9933   0.09151     2.750
NumLiverMetsGroup≥2        2.9583     0.3380   0.49889    17.543
SizeLiverMetsmmGroup≥50    5.8302     0.1715   0.38808    87.588
RASMutant                  4.9288     0.2029   0.66321    36.630
Oxaliplatin.HistoryYes         NA         NA        NA        NA

Concordance= 0.824  (se = 0.07 )
Likelihood ratio test= 10.15  on 6 df,   p=0.1
Wald test            = 7.56  on 6 df,   p=0.3
Score (logrank) test = 9.31  on 6 df,   p=0.2

#DFS by ACT treatment in MRD negative - Metachronous metastatic site

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Metachronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 76     28     NA    25.6      NA
ACT=TRUE  16      5     NA    14.6      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | Metachronous Met Site", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000      20.0000      26.0000       0.6142       0.0632       0.4784       0.7246 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        6.000        5.000        0.675        0.121        0.384        0.851 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 92, number of events= 33 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.2494    1.2833   0.4859 0.513    0.608

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     1.283     0.7792    0.4951     3.326

Concordance= 0.514  (se = 0.033 )
Likelihood ratio test= 0.28  on 1 df,   p=0.6
Wald test            = 0.26  on 1 df,   p=0.6
Score (logrank) test = 0.26  on 1 df,   p=0.6

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.28 (0.5-3.33); p = 0.608"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
                   1 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", dfs_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months 
                             1 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$Mets.Type=="Metachronous",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + NumLiverMetsGroup + 
    SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 92, number of events= 33 

                            coef exp(coef) se(coef)      z Pr(>|z|)  
ACTTRUE                 -0.21989   0.80261  0.50394 -0.436   0.6626  
GenderMale               0.83345   2.30124  0.47385  1.759   0.0786 .
Age.Group>70             0.49048   1.63310  0.37182  1.319   0.1871  
NumLiverMetsGroup≥2      0.05637   1.05799  0.42226  0.133   0.8938  
SizeLiverMetsmmGroup≥50  0.39198   1.47990  0.75395  0.520   0.6031  
RASMutant                0.23766   1.26827  0.39924  0.595   0.5517  
Oxaliplatin.HistoryYes   0.41547   1.51508  0.39143  1.061   0.2885  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.8026     1.2459    0.2989     2.155
GenderMale                 2.3012     0.4345    0.9091     5.825
Age.Group>70               1.6331     0.6123    0.7880     3.385
NumLiverMetsGroup≥2        1.0580     0.9452    0.4624     2.421
SizeLiverMetsmmGroup≥50    1.4799     0.6757    0.3377     6.486
RASMutant                  1.2683     0.7885    0.5799     2.774
Oxaliplatin.HistoryYes     1.5151     0.6600    0.7035     3.263

Concordance= 0.623  (se = 0.046 )
Likelihood ratio test= 6.09  on 7 df,   p=0.5
Wald test            = 5.71  on 7 df,   p=0.6
Score (logrank) test = 5.8  on 7 df,   p=0.6

#DFS by ACT treatment in MRD negative - 1 Liver Met

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="1",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 69     22     NA      31      NA
ACT=TRUE  22      6     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | 1 Liver Met", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000       21.000       20.000        0.674        0.063        0.534        0.780 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000       11.000        6.000        0.693        0.107        0.434        0.851 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 91, number of events= 28 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 0.3365    1.4000   0.4618 0.729    0.466

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE       1.4     0.7143    0.5663     3.461

Concordance= 0.536  (se = 0.038 )
Likelihood ratio test= 0.56  on 1 df,   p=0.5
Wald test            = 0.53  on 1 df,   p=0.5
Score (logrank) test = 0.54  on 1 df,   p=0.5

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 1.4 (0.57-3.46); p = 0.466"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
                   1 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", dfs_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months 
                             1 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="1",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + Mets.Type + SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + Mets.Type + 
    SizeLiverMetsmmGroup + RAS + Oxaliplatin.History, data = circ_data)

  n= 91, number of events= 28 

                            coef exp(coef) se(coef)      z Pr(>|z|)  
ACTTRUE                 -0.02522   0.97509  0.52282 -0.048   0.9615  
GenderMale               0.82189   2.27480  0.49890  1.647   0.0995 .
Age.Group>70             0.05878   1.06054  0.42162  0.139   0.8891  
Mets.TypeSynchronous    -0.38867   0.67796  0.55840 -0.696   0.4864  
SizeLiverMetsmmGroup≥50  1.26344   3.53757  0.63719  1.983   0.0474 *
RASMutant                0.82093   2.27261  0.45519  1.803   0.0713 .
Oxaliplatin.HistoryYes   0.25359   1.28864  0.50661  0.501   0.6167  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                        exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                    0.9751     1.0255    0.3500     2.717
GenderMale                 2.2748     0.4396    0.8556     6.048
Age.Group>70               1.0605     0.9429    0.4641     2.423
Mets.TypeSynchronous       0.6780     1.4750    0.2269     2.025
SizeLiverMetsmmGroup≥50    3.5376     0.2827    1.0147    12.333
RASMutant                  2.2726     0.4400    0.9312     5.546
Oxaliplatin.HistoryYes     1.2886     0.7760    0.4774     3.478

Concordance= 0.66  (se = 0.053 )
Likelihood ratio test= 8.34  on 7 df,   p=0.3
Wald test            = 9.02  on 7 df,   p=0.3
Score (logrank) test = 9.93  on 7 df,   p=0.2

#DFS by ACT treatment in MRD negative - 2≥ Liver Mets

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="≥2",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ACT, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ACT, data = circ_data)

           n events median 0.95LCL 0.95UCL
ACT=FALSE 23     11   20.3    6.15      NA
ACT=TRUE  11      2     NA      NA      NA

event_summary <- circ_data %>%
  group_by(ACT) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ACT, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("red","blue"), title="DFS - ctDNA MRD Negative ACT vs Observation | ≥2 Liver Met", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Observation", "ACT"), legend.title="")

summary(KM_curve, times= c(24))

Call: survfit(formula = surv_object ~ ACT, data = circ_data, conf.int = 0.95, 
    conf.type = "log-log")

                ACT=FALSE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        4.000       11.000        0.458        0.122        0.218        0.669 

                ACT=TRUE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
      24.000        4.000        2.000        0.818        0.116        0.447        0.951 

circ_data$ACT <- factor(circ_data$ACT, levels=c("TRUE","FALSE"))
cox_fit <- coxph(surv_object ~ ACT, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT, data = circ_data)

  n= 34, number of events= 13 

           coef exp(coef) se(coef)     z Pr(>|z|)
ACTFALSE 1.1484    3.1533   0.7692 1.493    0.135

         exp(coef) exp(-coef) lower .95 upper .95
ACTFALSE     3.153     0.3171    0.6983     14.24

Concordance= 0.602  (se = 0.062 )
Likelihood ratio test= 2.84  on 1 df,   p=0.09
Wald test            = 2.23  on 1 df,   p=0.1
Score (logrank) test = 2.48  on 1 df,   p=0.1

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 3.15 (0.7-14.24); p = 0.135"

#Fisher test for DFS percentages at 24
dfs_times <- c(24)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ACT == "TRUE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ACT == "FALSE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ACT == "TRUE")
  pos_total <- sum(circ_data$ACT == "FALSE")
  
  neg_surv <- neg_total - sum(circ_data$ACT == "TRUE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ACT == "FALSE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months 
            0.139822 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", dfs_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months 
                      0.139822 

#Adjusted HR "ACT vs no ACT" - age, gender, ECOG and pathological stage - ACT as reference
rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ctDNA.MRD=="NEGATIVE",]
circ_data <- circ_data[circ_data$NumLiverMetsGroup=="≥2",]
circ_data$DFS.months=circ_data$DFS.months-2
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ACT + Gender + Age.Group + Mets.Type + RAS + Oxaliplatin.History, data=circ_data) 
summary(cox_fit)

Call:
coxph(formula = surv_object ~ ACT + Gender + Age.Group + Mets.Type + 
    RAS + Oxaliplatin.History, data = circ_data)

  n= 34, number of events= 13 

                            coef exp(coef)  se(coef)      z Pr(>|z|)
ACTTRUE                -1.350297  0.259163  0.833079 -1.621    0.105
GenderMale             -0.588374  0.555229  0.792352 -0.743    0.458
Age.Group>70            0.399882  1.491649  0.597240  0.670    0.503
Mets.TypeSynchronous   -0.006411  0.993610  0.794551 -0.008    0.994
RASMutant              -0.396021  0.672993  0.629514 -0.629    0.529
Oxaliplatin.HistoryYes  0.149209  1.160916  0.690752  0.216    0.829

                       exp(coef) exp(-coef) lower .95 upper .95
ACTTRUE                   0.2592     3.8586   0.05064     1.326
GenderMale                0.5552     1.8011   0.11750     2.624
Age.Group>70              1.4916     0.6704   0.46269     4.809
Mets.TypeSynchronous      0.9936     1.0064   0.20936     4.716
RASMutant                 0.6730     1.4859   0.19596     2.311
Oxaliplatin.HistoryYes    1.1609     0.8614   0.29980     4.495

Concordance= 0.661  (se = 0.086 )
Likelihood ratio test= 4.56  on 6 df,   p=0.6
Wald test            = 3.86  on 6 df,   p=0.7
Score (logrank) test = 4.16  on 6 df,   p=0.7

#DFS by ctDNA Clearance ACT-treated at 3 months

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ACT==TRUE,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.MRD == "POSITIVE" & ctDNA.3months == "NEGATIVE" ~ 1,
    ctDNA.MRD == "POSITIVE" & ctDNA.3months == "POSITIVE" ~ 2
  ))

circ_data <- circ_data[circ_data$DFS.3mo.months>=0,]
survfit(Surv(time = circ_data$DFS.3mo.months, event = circ_data$DFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.3mo.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

   35 observations deleted due to missingness 
                 n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 8      4  14.00   13.80      NA
ctDNA.Dynamics=2 5      5   3.42    3.02      NA

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.3mo.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA Clearance from MRD to 3 months ACT-treated", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Clearance", "No Clearance"), legend.title="")

summary(KM_curve, times= c(6, 12, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

35 observations deleted due to missingness 
                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    6      7       1    0.875   0.117        0.387        0.981
   12      7       0    0.875   0.117        0.387        0.981
   24      3       3    0.500   0.177        0.152        0.775

                ctDNA.Dynamics=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
       6.000        2.000        3.000        0.400        0.219        0.052        0.753 

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 13, number of events= 9 
   (35 observations deleted due to missingness)

                             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.DynamicsNo Clearance  2.717    15.131    1.117 2.432    0.015 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance     15.13    0.06609     1.695     135.1

Concordance= 0.736  (se = 0.053 )
Likelihood ratio test= 8.61  on 1 df,   p=0.003
Wald test            = 5.92  on 1 df,   p=0.02
Score (logrank) test = 9.83  on 1 df,   p=0.002

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 15.13 (1.69-135.08); p = 0.015"

#Levels of MRD MTM/mL in Clearance at 3 months log10 transformation

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ACT==TRUE,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Clearance <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Clearance = case_when(
    ctDNA.MRD == "POSITIVE" & ctDNA.3months == "NEGATIVE" ~ 1,
    ctDNA.MRD == "POSITIVE" & ctDNA.3months == "POSITIVE" ~ 2
  ))

circ_data <- circ_data[circ_data$DFS.3mo.months>=0,]
circ_data <- subset(circ_data, !is.na(ctDNA.Clearance))

# Transform p_MRD_MTM with log10
circ_data$p_MRD_MTM <- as.numeric(as.character(circ_data$p_MRD_MTM))
circ_data$ctDNA.Clearance <- factor(circ_data$ctDNA.Clearance, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
median_p_MRD_MTM <- aggregate(p_MRD_MTM ~ ctDNA.Clearance, data = circ_data, FUN = median)
print(median_p_MRD_MTM)

# Create violin plot with log10 scale on y-axis
ggplot(circ_data, aes(x=ctDNA.Clearance, y=p_MRD_MTM, fill=ctDNA.Clearance)) +
  geom_violin(trim=FALSE) +
  scale_fill_manual(values=c("Clearance"="lightblue", "No Clearance"="salmon")) +
  geom_boxplot(width=0.1, fill="white", colour="black", alpha=0.5) +
  scale_y_log10(breaks=c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000)) +
  labs(title="MRD MTM/mL | Clearance at 3 months", x="Clearance at 3 months", y="MRD MTM/mL") +
  theme_minimal() +
  theme(legend.position="none")

m3_1v2 <- wilcox.test(p_MRD_MTM ~ ctDNA.Clearance,
                      data = circ_data[circ_data$ctDNA.Clearance %in% c("Clearance", "No Clearance"), ],
                      na.rm = TRUE)
print(m3_1v2)

    Wilcoxon rank sum exact test

data:  p_MRD_MTM by ctDNA.Clearance
W = 13, p-value = 0.3543
alternative hypothesis: true location shift is not equal to 0

#DFS by ctDNA Clearance ACT-treated at 6 months

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ACT==TRUE,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Dynamics <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Dynamics = case_when(
    ctDNA.MRD == "POSITIVE" & ctDNA.6months == "NEGATIVE" ~ 1,
    ctDNA.MRD == "POSITIVE" & ctDNA.6months == "POSITIVE" ~ 2
  ))

circ_data <- circ_data[circ_data$DFS.6mo.months>=0,]
survfit(Surv(time = circ_data$DFS.6mo.months, event = circ_data$DFS.Event)~ctDNA.Dynamics, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.6mo.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Dynamics, data = circ_data)

   36 observations deleted due to missingness 
                 n events median 0.95LCL 0.95UCL
ctDNA.Dynamics=1 9      5 11.696   9.429      NA
ctDNA.Dynamics=2 3      3  0.789   0.493      NA

event_summary <- circ_data %>%
  group_by(ctDNA.Dynamics) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.6mo.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Dynamics, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA Clearance from MRD to 6 months ACT-treated", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("Clearance", "No Clearance"), legend.title="")

summary(KM_curve, times= c(6, 24))

Call: survfit(formula = surv_object ~ ctDNA.Dynamics, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

36 observations deleted due to missingness 
                ctDNA.Dynamics=1 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    6      7       2    0.778   0.139        0.365        0.939
   24      2       3    0.444   0.166        0.136        0.719

                ctDNA.Dynamics=2 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     6.00000      1.00000      2.00000      0.33333      0.27217      0.00896      0.77415 

circ_data$ctDNA.Dynamics <- factor(circ_data$ctDNA.Dynamics, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
cox_fit <- coxph(surv_object ~ ctDNA.Dynamics, data=circ_data) 
ggforest(cox_fit,data = circ_data) 

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Dynamics, data = circ_data)

  n= 12, number of events= 8 
   (36 observations deleted due to missingness)

                             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.DynamicsNo Clearance 1.9221    6.8355   0.9345 2.057   0.0397 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.DynamicsNo Clearance     6.836     0.1463     1.095     42.68

Concordance= 0.658  (se = 0.073 )
Likelihood ratio test= 4.19  on 1 df,   p=0.04
Wald test            = 4.23  on 1 df,   p=0.04
Score (logrank) test = 5.53  on 1 df,   p=0.02

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 6.84 (1.09-42.68); p = 0.04"

#Levels of MRD MTM/mL in Clearance at 6 months log10 transformation

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$ACT==TRUE,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Clearance <- NA #first we create the variable for the ctDNA & NAC combination, and we assign values
circ_data <- circ_data %>%
  mutate(ctDNA.Clearance = case_when(
    ctDNA.MRD == "POSITIVE" & ctDNA.6months == "NEGATIVE" ~ 1,
    ctDNA.MRD == "POSITIVE" & ctDNA.6months == "POSITIVE" ~ 2
  ))

circ_data <- circ_data[circ_data$DFS.6mo.months>=0,]
circ_data <- subset(circ_data, !is.na(ctDNA.Clearance))

# Transform p_MRD_MTM with log10
circ_data$p_MRD_MTM <- as.numeric(as.character(circ_data$p_MRD_MTM))
circ_data$ctDNA.Clearance <- factor(circ_data$ctDNA.Clearance, levels=c("1","2"), labels = c("Clearance", "No Clearance"))
median_p_MRD_MTM <- aggregate(p_MRD_MTM ~ ctDNA.Clearance, data = circ_data, FUN = median)
print(median_p_MRD_MTM)

# Create violin plot with log10 scale on y-axis
ggplot(circ_data, aes(x=ctDNA.Clearance, y=p_MRD_MTM, fill=ctDNA.Clearance)) +
  geom_violin(trim=FALSE) +
  scale_fill_manual(values=c("Clearance"="lightblue", "No Clearance"="salmon")) +
  geom_boxplot(width=0.1, fill="white", colour="black", alpha=0.5) +
  scale_y_log10(breaks=c(0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000)) +
  labs(title="MRD MTM/mL | Clearance at 6 months", x="Clearance at 6 months", y="MRD MTM/mL") +
  theme_minimal() +
  theme(legend.position="none")

m3_1v2 <- wilcox.test(p_MRD_MTM ~ ctDNA.Clearance,
                      data = circ_data[circ_data$ctDNA.Clearance %in% c("Clearance", "No Clearance"), ],
                      na.rm = TRUE)
print(m3_1v2)

    Wilcoxon rank sum exact test

data:  p_MRD_MTM by ctDNA.Clearance
W = 16, p-value = 0.7273
alternative hypothesis: true location shift is not equal to 0

#DFS by ctDNA at the Surveillance Window - Landmark 10 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Surveillance!="",]
circ_data$DFS.months=circ_data$DFS.months-2.5
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)~ctDNA.Surveillance, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event) ~ 
    ctDNA.Surveillance, data = circ_data)

                             n events median 0.95LCL 0.95UCL
ctDNA.Surveillance=NEGATIVE 99     18     NA      NA      NA
ctDNA.Surveillance=POSITIVE 53     44   5.81    3.05    12.1

event_summary <- circ_data %>%
  group_by(ctDNA.Surveillance) %>%
  summarise(
    Total = n(),
    Events = sum(DFS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$DFS.months, event = circ_data$DFS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="DFS - ctDNA Surveillance window", ylab= "Disease-Free Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Surveillance=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     41      16    0.815  0.0434        0.711        0.884
   30     31       1    0.795  0.0468        0.684        0.870
   36     15       1    0.768  0.0522        0.646        0.853

                ctDNA.Surveillance=POSITIVE 
        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI 
     24.0000       2.0000      44.0000       0.0859       0.0501       0.0200       0.2142 

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance, data = circ_data)

  n= 152, number of events= 62 

                              coef exp(coef) se(coef)     z Pr(>|z|)    
ctDNA.SurveillancePOSITIVE  2.3644   10.6382   0.2982 7.929 2.21e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.SurveillancePOSITIVE     10.64      0.094      5.93     19.09

Concordance= 0.757  (se = 0.025 )
Likelihood ratio test= 73.09  on 1 df,   p=<2e-16
Wald test            = 62.87  on 1 df,   p=2e-15
Score (logrank) test = 91.54  on 1 df,   p=<2e-16

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 10.64 (5.93-19.09); p = 0"

#Fisher test for DFS percentages at 24, 30, and 36 months
dfs_times <- c(24, 30, 36)
p_values <- sapply(dfs_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.Surveillance == "NEGATIVE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.Surveillance == "POSITIVE" & circ_data$DFS.months >= time & circ_data$DFS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.Surveillance == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.Surveillance == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.Surveillance == "NEGATIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.Surveillance == "POSITIVE" & circ_data$DFS.Event == 1 & circ_data$DFS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", dfs_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
        4.109654e-16         1.518274e-15         4.230750e-15 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", dfs_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months Corrected p-value at 30 months Corrected p-value at 36 months 
                  1.232896e-15                   4.554823e-15                   1.269225e-14 

#OS by ctDNA at the Surveillance Window - Landmark 10 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Surveillance!="",]
circ_data$OS.months=circ_data$OS.months-2.5
circ_data <- circ_data[circ_data$OS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.months, event = circ_data$OS.Event)~ctDNA.Surveillance, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.months, event = circ_data$OS.Event) ~ 
    ctDNA.Surveillance, data = circ_data)

                             n events median 0.95LCL 0.95UCL
ctDNA.Surveillance=NEGATIVE 99      3     NA      NA      NA
ctDNA.Surveillance=POSITIVE 55      6   41.8    37.3      NA

event_summary <- circ_data %>%
  group_by(ctDNA.Surveillance) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.Surveillance, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - ctDNA Surveillance window", ylab= "Overall Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 30, 36))

Call: survfit(formula = surv_object ~ ctDNA.Surveillance, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.Surveillance=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     55       1    0.990  0.0104        0.928        0.999
   30     39       1    0.965  0.0259        0.855        0.992
   36     18       1    0.941  0.0351        0.817        0.982

                ctDNA.Surveillance=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     22       2    0.955  0.0319        0.827        0.989
   30     11       2    0.858  0.0715        0.641        0.948
   36      5       0    0.858  0.0715        0.641        0.948

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.Surveillance, data = circ_data)

  n= 154, number of events= 9 

                             coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.SurveillancePOSITIVE 1.6394    5.1523   0.7188 2.281   0.0226 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                           exp(coef) exp(-coef) lower .95 upper .95
ctDNA.SurveillancePOSITIVE     5.152     0.1941     1.259     21.08

Concordance= 0.678  (se = 0.094 )
Likelihood ratio test= 5.59  on 1 df,   p=0.02
Wald test            = 5.2  on 1 df,   p=0.02
Score (logrank) test = 6.33  on 1 df,   p=0.01

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 5.15 (1.26-21.08); p = 0.023"

#Fisher test for OS percentages at 24, 30, and 36 months
os_times <- c(24, 30, 36)
p_values <- sapply(os_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.months >= time & circ_data$OS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.months >= time & circ_data$OS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.Event == 1 & circ_data$OS.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.Event == 1 & circ_data$OS.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", os_times, " months")
print(p_values)

p-value at 24 months p-value at 30 months p-value at 36 months 
          0.12280702           0.02179403           0.04276279 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", os_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months Corrected p-value at 30 months Corrected p-value at 36 months 
                    0.36842105                     0.06538209                     0.12828838 

#Multivariate cox regression at Surveillance Window for DFS - Landmark 10 weeks

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.Surveillance!="",]
circ_data$DFS.months=circ_data$DFS.months-2.5
circ_data <- circ_data[circ_data$DFS.months>=0,]
circ_datadf <- as.data.frame(circ_data)

circ_data$ctDNA.Surveillance <- factor(circ_data$ctDNA.Surveillance, levels=c("NEGATIVE","POSITIVE"), labels = c("Negative", "Positive"))
circ_data$CEA.Baseline <- factor(circ_data$CEA.Baseline, levels=c("NEGATIVE","POSITIVE"), labels = c("Normal", "Elevated"))
circ_data$Gender <- factor(circ_data$Gender, levels = c("Female", "Male"))
circ_data$Age.Group <- factor(circ_data$Age.Group, levels = c("1", "2"), labels = c("<70", ">70"))
circ_data$PrimSite <- factor(circ_data$PrimSite, levels = c("Right-sided colon", "Left-sided colon"), labels = c("Right-sided", "Left-sided"))
circ_data$NumLiverMetsGroup <- factor(circ_data$NumLiverMetsGroup, levels = c("1", "≥2"))
circ_data$SizeLiverMetsmmGroup <- factor(circ_data$SizeLiverMetsmmGroup, levels = c("<50", "≥50"))
circ_data$ECOG <- factor(circ_data$ECOG, levels = c("0", "1"))
circ_data$pT <- factor(circ_data$pT, levels = c("T1-T2", "T3-T4"))
circ_data$pN <- factor(circ_data$pN, levels = c("N0", "N1-N2"))
circ_data$ACT <- factor(circ_data$ACT, levels = c("FALSE", "TRUE"), labels = c("Observation", "Chemotherapy"))
circ_data$Postop.Complication <- factor(circ_data$Postop.Complication, levels = c("No", "Yes"))
circ_data$Oxaliplatin.History <- factor(circ_data$Oxaliplatin.History, levels = c("No", "Yes"))
circ_data$Mets.Type <- factor(circ_data$Mets.Type, levels = c("Metachronous", "Synchronous"))
circ_data$MSI <- factor(circ_data$MSI, levels = c("MSS", "MSI-High"), labels = c("MSS", "MSI-High"))
circ_data$BRAF.V600E <- factor(circ_data$BRAF.V600E, levels = c("WT", "MUT"), labels = c("Wild-Type", "V600E"))
circ_data$RAS <- factor(circ_data$RAS, levels = c("WT", "MUT"), labels = c("Wild-Type", "Mutant"))
surv_object <- Surv(time = circ_data$DFS.MRD.months, event = circ_data$DFS.Event) 
cox_fit <- coxph(surv_object ~ ctDNA.Surveillance + CEA.Baseline + Age.Group + PrimSite + NumLiverMetsGroup + SizeLiverMetsmmGroup + pT + pN + RAS + Mets.Type + Postop.Complication + Oxaliplatin.History, data=circ_data) 
ggforest(cox_fit, data = circ_data, main = "Multivariate Regression Model for DFS", refLabel = "Reference Group")

test.ph <- cox.zph(cox_fit)

#OS by ctDNA at the MRD Window - pts with Radiological Recurrence

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]
circ_data <- circ_data[circ_data$ctDNA.MRD!="",]
circ_data <- circ_data[circ_data$OS.MRD.months>=0,]
circ_datadf <- as.data.frame(circ_data)

survfit(Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)~ctDNA.MRD, data = circ_data)

Call: survfit(formula = Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event) ~ 
    ctDNA.MRD, data = circ_data)

                    n events median 0.95LCL 0.95UCL
ctDNA.MRD=NEGATIVE 42      2     NA      NA      NA
ctDNA.MRD=POSITIVE 53     10   43.4      NA      NA

event_summary <- circ_data %>%
  group_by(ctDNA.MRD) %>%
  summarise(
    Total = n(),
    Events = sum(OS.Event),
    Fraction = Events / n(),
    Percentage = (Events / n()) * 100
  )
print(event_summary)

surv_object <-Surv(time = circ_data$OS.MRD.months, event = circ_data$OS.Event)
KM_curve <- survfit(surv_object ~ ctDNA.MRD, data = circ_data,conf.int=0.95,conf.type="log-log") 
ggsurvplot(KM_curve, data = circ_data, pval = FALSE, conf.int = FALSE, risk.table = TRUE, break.time.by=6, palette=c("blue","red"), title="OS - Radiological Recurrence | ctDNA MRD window", ylab= "Overall Survival", xlab="Time from Landmark Time point (Months)", legend.labs=c("ctDNA Negative", "ctDNA Positive"), legend.title="")

summary(KM_curve, times= c(24, 36))

Call: survfit(formula = surv_object ~ ctDNA.MRD, data = circ_data, 
    conf.int = 0.95, conf.type = "log-log")

                ctDNA.MRD=NEGATIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     25       0    1.000  0.0000           NA           NA
   36      8       1    0.933  0.0644        0.613         0.99

                ctDNA.MRD=POSITIVE 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
   24     18       6     0.82  0.0709        0.627        0.919
   36      4       3     0.68  0.0942        0.459        0.827

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels=c("NEGATIVE","POSITIVE"))
cox_fit <- coxph(surv_object ~ ctDNA.MRD, data=circ_data) 
ggforest(cox_fit,data = circ_data)

summary(cox_fit)

Call:
coxph(formula = surv_object ~ ctDNA.MRD, data = circ_data)

  n= 95, number of events= 12 

                    coef exp(coef) se(coef)     z Pr(>|z|)  
ctDNA.MRDPOSITIVE 1.8047    6.0781   0.7937 2.274    0.023 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                  exp(coef) exp(-coef) lower .95 upper .95
ctDNA.MRDPOSITIVE     6.078     0.1645     1.283      28.8

Concordance= 0.735  (se = 0.031 )
Likelihood ratio test= 6.88  on 1 df,   p=0.009
Wald test            = 5.17  on 1 df,   p=0.02
Score (logrank) test = 6.56  on 1 df,   p=0.01

cox_fit_summary <- summary(cox_fit)

# Extract values for HR, 95% CI, and p-value
HR <- cox_fit_summary$coefficients[2]
lower_CI <- cox_fit_summary$conf.int[3]
upper_CI <- cox_fit_summary$conf.int[4]
p_value <- cox_fit_summary$coefficients[5]
label_text <- paste0("HR = ", round(HR, 2), " (", round(lower_CI, 2), "-", round(upper_CI, 2), "); p = ", round(p_value, 3))
print(label_text)

[1] "HR = 6.08 (1.28-28.8); p = 0.023"

#Fisher test for OS percentages at 24 and 36 months
os_times <- c(24, 36)
p_values <- sapply(os_times, function(time) {
  neg_count <- sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  pos_count <- sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.MRD.months >= time & circ_data$OS.Event == 0)
  neg_total <- sum(circ_data$ctDNA.MRD == "NEGATIVE")
  pos_total <- sum(circ_data$ctDNA.MRD == "POSITIVE")
  
  neg_surv <- neg_total - sum(circ_data$ctDNA.MRD == "NEGATIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  pos_surv <- pos_total - sum(circ_data$ctDNA.MRD == "POSITIVE" & circ_data$OS.Event == 1 & circ_data$OS.MRD.months < time)
  
  surv_matrix <- matrix(c(neg_surv, pos_surv, neg_total - neg_surv, pos_total - pos_surv), nrow = 2)
  test_result <- fisher.test(surv_matrix)
  return(test_result$p.value)
})
names(p_values) <- paste0("p-value at ", os_times, " months")
print(p_values)

p-value at 24 months p-value at 36 months 
          0.03245083           0.03892912 

p_values_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_values_adjusted) <- paste0("Corrected p-value at ", os_times, " months")
print(p_values_adjusted)

Corrected p-value at 24 months Corrected p-value at 36 months 
                    0.06490165                     0.07785824 

#Percentage of ctDNA MRD Window positivity in pts with liver metastasis

rm(list = ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data %>%
  filter(LiverMets == "TRUE" & RFS.Event == "TRUE" & ctDNA.MRD != "")
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
positive_rate <- sum(circ_data$ctDNA.MRD == "Positive" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Positive") * 100
positive_ci <- binconf(sum(circ_data$ctDNA.MRD == "Positive" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Positive"),
                       alpha = 0.05)[c(2, 3)] * 100
negative_rate <- sum(circ_data$ctDNA.MRD == "Negative" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Negative") * 100
negative_ci <- binconf(sum(circ_data$ctDNA.MRD == "Negative" & grepl("liver", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Negative"),
                       alpha = 0.05)[c(2, 3)] * 100
data <- data.frame(
  ctDNA.MRD = c("Positive", "Negative"),
  percentage = c(positive_rate, negative_rate),
  lower_ci = c(positive_ci[1], negative_ci[1]),
  upper_ci = c(positive_ci[2], negative_ci[2])
)
cross_tabs <- list(
  table(circ_data$ctDNA.MRD, grepl("liver", circ_data$RelSite, ignore.case = TRUE)),
  table(circ_data$ctDNA.MRD, grepl("lung", circ_data$RelSite, ignore.case = TRUE))  # Example for another site
)
chi_tests <- lapply(cross_tabs, chisq.test)
p_values <- sapply(chi_tests, function(test) test$p.value)
num_tests <- length(p_values)  # Number of tests to be adjusted for
p_values_adjusted <- p.adjust(p_values, method = "bonferroni", n = num_tests)
print(data)

print(cross_tabs)

[[1]]
          
           FALSE TRUE
  Negative    28   14
  Positive    14   39

[[2]]
          
           FALSE TRUE
  Negative    16   26
  Positive    39   14

print(chi_tests)

[[1]]

    Pearson's Chi-squared test with Yates' continuity correction

data:  X[[i]]
X-squared = 13.803, df = 1, p-value = 0.000203


[[2]]

    Pearson's Chi-squared test with Yates' continuity correction

data:  X[[i]]
X-squared = 10.695, df = 1, p-value = 0.001074

print(paste("Original p-values:", format.pval(p_values, digits = 4)))

[1] "Original p-values: 0.0002" "Original p-values: 0.0011"

print(paste("Adjusted p-values (Bonferroni):", format.pval(p_values_adjusted, digits = 4)))

[1] "Adjusted p-values (Bonferroni): 0.0004" "Adjusted p-values (Bonferroni): 0.0021"

barplot <- ggplot(data, aes(x = ctDNA.MRD, y = percentage, fill = ctDNA.MRD)) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2) +
  geom_text(aes(label = paste0(round(percentage, 1), "%")), vjust = -0.5) +
  labs(
    x = "ctDNA status at the MRD status",
    y = "Proportion of patients with recurrence in the liver",
    caption = paste("Chi-squared test p-value: ", format.pval(p_values[1], digits = 4), 
                    "\nAdjusted p-value (Bonferroni): ", format.pval(p_values_adjusted[1], digits = 4))
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 100)) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) +
  theme_minimal()
print(barplot)

#Percentage of ctDNA MRD Window positivity in pts with lung metastasis

rm(list = ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data %>%
  filter(LiverMets == "TRUE" & RFS.Event == "TRUE" & ctDNA.MRD != "")
circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
positive_rate <- sum(circ_data$ctDNA.MRD == "Positive" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Positive") * 100
positive_ci <- binconf(sum(circ_data$ctDNA.MRD == "Positive" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Positive"),
                       alpha = 0.05)[c(2, 3)] * 100
negative_rate <- sum(circ_data$ctDNA.MRD == "Negative" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)) / sum(circ_data$ctDNA.MRD == "Negative") * 100
negative_ci <- binconf(sum(circ_data$ctDNA.MRD == "Negative" & grepl("lung", circ_data$RelSite, ignore.case = TRUE)),
                       sum(circ_data$ctDNA.MRD == "Negative"),
                       alpha = 0.05)[c(2, 3)] * 100
data <- data.frame(
  ctDNA.MRD = c("Positive", "Negative"),
  percentage = c(positive_rate, negative_rate),
  lower_ci = c(positive_ci[1], negative_ci[1]),
  upper_ci = c(positive_ci[2], negative_ci[2])
)

cross_tabs <- list(
  table(circ_data$ctDNA.MRD, grepl("lung", circ_data$RelSite, ignore.case = TRUE)),
  table(circ_data$ctDNA.MRD, grepl("liver", circ_data$RelSite, ignore.case = TRUE))
)
chi_tests <- lapply(cross_tabs, chisq.test)
p_values <- sapply(chi_tests, function(test) test$p.value)
num_tests <- length(p_values)  # Number of tests to be adjusted for
p_values_adjusted <- p.adjust(p_values, method = "bonferroni", n = num_tests)
print(data)

print(cross_tabs)

[[1]]
          
           FALSE TRUE
  Negative    16   26
  Positive    39   14

[[2]]
          
           FALSE TRUE
  Negative    28   14
  Positive    14   39

print(chi_tests)

[[1]]

    Pearson's Chi-squared test with Yates' continuity correction

data:  X[[i]]
X-squared = 10.695, df = 1, p-value = 0.001074


[[2]]

    Pearson's Chi-squared test with Yates' continuity correction

data:  X[[i]]
X-squared = 13.803, df = 1, p-value = 0.000203

print(paste("Original p-values:", format.pval(p_values, digits = 4)))

[1] "Original p-values: 0.0011" "Original p-values: 0.0002"

print(paste("Adjusted p-values (Bonferroni):", format.pval(p_values_adjusted, digits = 4)))

[1] "Adjusted p-values (Bonferroni): 0.0021" "Adjusted p-values (Bonferroni): 0.0004"

barplot <- ggplot(data, aes(x = ctDNA.MRD, y = percentage, fill = ctDNA.MRD)) +
  geom_bar(stat = "identity") +
  geom_errorbar(aes(ymin = lower_ci, ymax = upper_ci), width = 0.2) +
  geom_text(aes(label = paste0(round(percentage, 1), "%")), vjust = -0.5) +
  labs(
    x = "ctDNA status at the MRD status",
    y = "Proportion of patients with recurrence in the lung",
    caption = paste("Chi-squared test p-value: ", format.pval(p_values[1], digits = 4), 
                    "\nAdjusted p-value (Bonferroni): ", format.pval(p_values_adjusted[1], digits = 4))
  ) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 100)) +
  scale_fill_manual(values = c("Negative" = "blue", "Positive" = "red")) +
  theme_minimal()
print(barplot)

#Barplot with Recurrence Sites (Liver vs Others) by ctDNA at the MRD Window

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]

circ_data$ctDNA.MRD <- factor(circ_data$ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE"), labels = c("Negative", "Positive"))
circ_data$Liver <- factor(circ_data$Liver, levels = c("FALSE", "TRUE"), labels = c("Others", "Liver"))
contingency_table <- table(circ_data$ctDNA.MRD, circ_data$Liver)
chi_square_test <- chisq.test(contingency_table)
print(chi_square_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  contingency_table
X-squared = 18.352, df = 1, p-value = 1.836e-05

fisher_exact_test <- fisher.test(contingency_table)
print(fisher_exact_test)

    Fisher's Exact Test for Count Data

data:  contingency_table
p-value = 9.391e-06
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  2.744866 20.863460
sample estimates:
odds ratio 
  7.311049 

print(contingency_table)

          
           Others Liver
  Negative     30    12
  Positive     13    39

p_values <- c(chi_square_test$p.value, fisher_exact_test$p.value)
p_adjusted <- p.adjust(p_values, method = "bonferroni")
names(p_adjusted) <- c("Chi-Square Test", "Fisher's Exact Test")
print(p_adjusted)

    Chi-Square Test Fisher's Exact Test 
       3.672023e-05        1.878209e-05 

table_df <- as.data.frame(contingency_table)
table_df$Total <- ave(table_df$Freq, table_df$Var1, FUN = sum)
table_df$Percentage <- table_df$Freq / table_df$Total
table_df$MiddlePercentage <- table_df$Percentage / 2
ggplot(table_df, aes(x = Var1, y = Percentage, fill = Var2)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = MiddlePercentage, label = Freq), position = "stack", color = "black", vjust = 1.5, size = 7) +
  theme_minimal() +
  labs(title = "Patients with Radiological Recurrence", 
       x = "ctDNA at the MRD Window", 
       y = "Patients (%)", 
       fill = "Recurrence Site",
       caption = paste("Chi-squared test p-value: ", format.pval(chi_square_test$p.value), "; Bonferroni corrected p-value: ", format.pval(p_adjusted["Chi-Square Test"]))) +
  scale_y_continuous(labels = scales::percent_format()) +
  scale_fill_manual(values = c("Others" = "blue", "Liver" = "red")) + # define custom colors
  theme(axis.text.x = element_text(angle = 0, hjust = 1.5, size = 14), # increase x-axis text size
        axis.text.y = element_text(size = 14, color = "black"), # increase y-axis text size
        axis.title.x = element_text(size = 14, color = "black"), # increase x-axis label size
        axis.title.y = element_text(size = 14, color = "black"), # increase y-axis label size
        legend.text = element_text(size = 12, color = "black"))  # increase Recurrence label size

#Table with recurrence sites by ctDNA at the MRD Window

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]

circ_data_subset1 <- circ_data %>%
  select(
    Rec.Site) %>%
  mutate(
    Rec.Site = factor(Rec.Site, levels = c("Local/LN", "Liver", "Lung", "Peritoneum", "Peritoneum & Others", "Brain")))

circ_data1 <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data1 <- circ_data1[circ_data1$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]

circ_data_subset2 <- circ_data %>%
  select(
    Rec.Site,
    ctDNA.MRD) %>%
  mutate(
    Rec.Site = factor(Rec.Site, levels = c("Local/LN", "Liver", "Lung", "Peritoneum", "Peritoneum & Others", "Brain")),
    ctDNA.MRD = factor(ctDNA.MRD, levels = c("NEGATIVE", "POSITIVE")))
Overall <- circ_data_subset1 %>%
  tbl_summary(
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  bold_labels()
Overall

Characteristic	N = 951
Rec.Site
Local/LN	4 (4.3%)
Liver	48 (51%)
Lung	30 (32%)
Peritoneum	4 (4.3%)
Peritoneum & Others	7 (7.4%)
Brain	1 (1.1%)
Unknown	1
1 n (%)

ByctDNA_MRD <- circ_data_subset2 %>%
  tbl_summary(
    by = ctDNA.MRD, # Subgroup by ctDNA.MRD
    statistic = list(
      all_continuous() ~ "{median} ({min} - {max})",
      all_categorical() ~ "{n} ({p}%)")) %>%
  add_p(test = all_categorical() ~ "fisher.test", pvalue_fun = ~style_pvalue(p.adjust(.x, method = "bonferroni"))) %>%
  bold_labels()
ByctDNA_MRD

Characteristic	NEGATIVE, N = 421	POSITIVE, N = 531	p-value2
Rec.Site			<0.001
Local/LN	1 (2.4%)	3 (5.8%)
Liver	12 (29%)	36 (69%)
Lung	23 (55%)	7 (13%)
Peritoneum	3 (7.1%)	1 (1.9%)
Peritoneum & Others	3 (7.1%)	4 (7.7%)
Brain	0 (0%)	1 (1.9%)
Unknown	0	1
1 n (%)
2 Fisher’s exact test

merged_table <- tbl_merge(tbls=list(Overall, ByctDNA_MRD))
merged_table

Characteristic	Table 1	Table 2
	N = 951	NEGATIVE, N = 421	POSITIVE, N = 531	p-value2
Rec.Site				<0.001
Local/LN	4 (4.3%)	1 (2.4%)	3 (5.8%)
Liver	48 (51%)	12 (29%)	36 (69%)
Lung	30 (32%)	23 (55%)	7 (13%)
Peritoneum	4 (4.3%)	3 (7.1%)	1 (1.9%)
Peritoneum & Others	7 (7.4%)	3 (7.1%)	4 (7.7%)
Brain	1 (1.1%)	0 (0%)	1 (1.9%)
Unknown	1	0	1
1 n (%)
2 Fisher’s exact test

fit1 <- as_flex_table(
  merged_table,
  include = everything(),
  return_calls = FALSE,
  strip_md_bold = TRUE)
fit1

	Table 1	Table 2
Characteristic	N = 951	NEGATIVE, N = 421	POSITIVE, N = 531	p-value2
Rec.Site				<0.001
Local/LN	4 (4.3%)	1 (2.4%)	3 (5.8%)
Liver	48 (51%)	12 (29%)	36 (69%)
Lung	30 (32%)	23 (55%)	7 (13%)
Peritoneum	4 (4.3%)	3 (7.1%)	1 (1.9%)
Peritoneum & Others	7 (7.4%)	3 (7.1%)	4 (7.7%)
Brain	1 (1.1%)	0 (0%)	1 (1.9%)
Unknown	1	0	1
1n (%)
2Fisher's exact test

save_as_docx(fit1, path= "~/Downloads/merged_table.docx")

#Detection ctDNA rates based on sites of relapse

rm(list=ls())
setwd("~/Downloads")
circ_data <- read.csv("Galaxy Liver Mets Data_20240626 Complete Dataset.csv")
circ_data <- circ_data[circ_data$LiverMets=="TRUE",]
circ_data <- circ_data[circ_data$RFS.Event=="TRUE",]
circ_data <- subset(circ_data, !is.na(Rec.Site))

# Create a table of counts for the "Rec.Site" variable
relsite_counts <- table(circ_data$Rec.Site)
relsite_df <- as.data.frame(relsite_counts)
names(relsite_df) <- c("RelSite", "Count")
circ_data_pos_mrd <- circ_data[circ_data$ctDNA.MRD=="POSITIVE",]
circ_data_pos_anytime <- circ_data[circ_data$ctDNA.anytime=="POSITIVE",]
pos_counts_mrd <- table(circ_data_pos_mrd$Rec.Site)
pos_counts_anytime <- table(circ_data_pos_anytime$Rec.Site)
relsite_df$MRDPos_Count <- ifelse(is.na(match(relsite_df$RelSite, names(pos_counts_mrd))), 0, pos_counts_mrd[match(relsite_df$RelSite, names(pos_counts_mrd))])
relsite_df$MRDPos_Count[is.na(relsite_df$MRDPos_Count)] <- 0
relsite_df$AnytimePos_Count <- ifelse(is.na(match(relsite_df$RelSite, names(pos_counts_anytime))), 0, pos_counts_anytime[match(relsite_df$RelSite, names(pos_counts_anytime))])
relsite_df$AnytimePos_Count[is.na(relsite_df$AnytimePos_Count)] <- 0
relsite_df$Percent <- (relsite_df$Count / sum(relsite_df$Count)) * 100
relsite_df$MRDPos_Percent <- (relsite_df$MRDPos_Count / relsite_df$Count) * 100
relsite_df$AnytimePos_Percent <- (relsite_df$AnytimePos_Count / relsite_df$Count) * 100
total_observations <- sum(relsite_df$Count)
total_pos_mrd <- sum(relsite_df$MRDPos_Count)
total_pos_anytime <- sum(relsite_df$AnytimePos_Count)
total_row <- data.frame(RelSite = "Total", Count = total_observations, MRDPos_Count = total_pos_mrd, AnytimePos_Count = total_pos_anytime, Percent = 100, MRDPos_Percent = (total_pos_mrd / total_observations) * 100, AnytimePos_Percent = (total_pos_anytime / total_observations) * 100)
relsite_df <- rbind(relsite_df, total_row)
print(relsite_df)

ft <- flextable(relsite_df)
doc <- read_docx() %>%
  body_add_flextable(value = ft)
print(doc, target = "relsite_df.docx")